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Abstract 



Relaxation phenomena in glasses can be related to jump processes between 
different minima of the potential energy in the configuration space. These 
transitions play a key role in the low temperature regime, giving rise to tun- 
neling systems responsible for the anomalous specific heat and thermal con- 
ductivity in disordered solids with respect to crystals. By using a recently 
developed numerical algorithm, we study the potential energy landscape of 
silica clusters, taking as a starting point the location of first order saddle 
points. This allows us to find a great number of adjacent minima. We an- 
alyze the degree of cooperativity of these transitions and the connection of 
physical properties with the topology of the configuration space. We also 
identify two-level systems (pairs of minima constituting a tunneling system) 
and calculate the quantum mechanical ground state splitting by means of the 
WKB approximation. 

PACS numbers : 61.43.Fs, 64.70.Pf, 82.20.Wt 
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I. INTRODUCTION 



The theoretical investigation of the properties of disordered systems is a very difficult 
task, because the lack of simmetry usually prevents analytical approaches. Many authors 
have found it convenient to rely on the concept of potential energy landscape in the configu- 
ration space (see e.g. Frauenfelder et al. 1991, Berry 1993, Heuer and Silbcy 1993, Mohanty 
et al. 1994, Heuer 1997, Angelani et al. 1998, Mousseau and Barkema 1998, Wales et al 
1998). 

Prom a purely qualitative point of view, we may imagine the multidimensional analogue 
of a surface rich in climbing points (first order saddles) and valleys (the various minima) . A 
detailed analysis of the energy landscape topology requires numerical simulation. 

In this work we analyze the potential energy landscape of Si02 clusters, namely [SiO2]20; 
[Si02]3o, and [Si02]5o, which are comparable in size to (or larger than) monoatomic and 
binary systems so far investigated. After a short description of the numerical methods used 
to move up and down the multidimensional hypersurface, we report our results regarding: 

1. the structure of the clusters (resembling that of the bulk solid, although affected by 
surface effects); 

2. the energy distribution of the stationary points of the potential energy function; 

3. a physical interpretation of geometrical quantities uselful to characterize transitions 
between different minima; 

4. the identification of tunneling systems. 



II. NUMERICAL METHODS 

We simulate the interaction among Si and O atoms with the recently developed pair 
potential by Van Beest et al. (1990), modified with a very short range contribution (Guissani 
and Guillot 1996), necessary to avoid unphysical divergencies. It is given by 



^ij = — + ttij exp {-hifij) - ^ + 4eij 



(Tij 



24 



(1) 



where both i and j indexes run on all Si and O atoms, and the values of the parameters 
were determined by Van Beest et al. (1990) and Guissani and Guillot (1996). This potential 
has already been widely used (with slight modifications) in molecular dynamics studies of 
the liquid and glassy phases of silica (sec e.g. VoUmair et al. 1996, Horbach et al. 1996, 
Taraskin and Elliott 1997). As our purpose is to simulate clusters, we used free boundary 
conditions. We adopted the procedure described in detail in Daldoss et al. (1998, 1999) 
to locate minimum-saddle-minimum triplets, that are the key ingredients to describe jump 
processes. Schematically the steps are: 

• descent towards a minimum by the conjugate gradient method (Press et al. 1986), 
starting from a randomly chosen configuration; 
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• ascent towards the vicinity of a saddle, following the eigenvector corresponding to the 
lowest eigenvalue; 

• once the potential energy along the path of the previous item starts decreasing, we 
take the corresponding configuration {hill- climbing point, Berry 1993) as the starting 
point for a Newton-Raphson search (Press et al. 1986) of the first order saddle point; 

• a descent on both sides of the saddle following again the eigenvector of the minimum 
eigenvalue. 

This technique provides approximate adjacent minima that are subsequently fed into the 
Newton algorithm for accurate location. With this numerical procedure we found many 
thousand minimum-saddle- minimum triplets (hereafter double well potentials, DWP). 

III. RESULTS 

A. The Structure of the clusters 

We obtain information on the structure of the cluster by the radial distribution function 
g{r), defined in the same way as for extended systems: it shows the bond lengths, the short 
range order, and also allows for the calculation of the co-ordination number. The results are 
reported in Fig.l, which was obtained by averaging over the various configurations obtained 
(both minima and saddles) to have good statistics. We have plotted the partial g{r) for the 
different bonds (Si-Si, Si-0, 0-0) and cluster sizes. The resulting bond lengths obtained by 
the main peaks are in good agreement with the experimental data on bulk structures (see 
table III) and with simulations with periodic boundary conditions (Taraskin and Elliott 1997). 
Nevertheless we observe smaller peaks on the low-distance side of the nearest-neighbour peak 
in the Si-Si and 0-0 bonds. In our opinion this is due to surface effects: in fact it appears 
from the figure that these anomalies tend to disappear with increasing system size. Very 
satisfactory is also the result concerning the co-ordination number: it is 4 for Si-0 bond, 
in agreement with the tetrahedral structure typical of Si02. This is also confirmed by the 
analysis of Si-Si, also equal to 4, and 0-0, equal to 6. 

B. Energy distributions of stationary points 

We now analyze the topological features of the potential energy hypersurface; at first 
we present the energy distribution of the stationary points (saddles and minima) that form 
double wells. In Fig. 2 we report the three energy distributions for saddles, lower minima and 
upper minima of the various DWP, respectively; only the result for 150 atoms is presented, 
the situation being quite analogous in the other two cases. We note that the distributions 
are very similar, with defined peaks superimposed onto a broad backgroud; the peaks are 
due to the fact that certain configurations are favoured for given values of the binding 
energy, and so they act like attraction basins during the descent towards the minima. The 
presence of these structures in the distribution is more evident with increasing system size: 
in [Si02]2o the curve looks smoother (Brangian 1998). It should be remarked that Fig. 2 
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does not refer to the total distribution of stationary points, because their number increases 
exponentially with the system size: we sampled partially the configuration space, and we 
are quite confident not to have introduced systematic errors in this sampling. The only 
thing we can note is that maybe our numerical investigation is in some way prevented from 
reaching very low lying (that is crystalline-like) configurations: since we are not interested in 
a careful thermodynamical analysis (Doye and Wales 1998), but only in transitions between 
stable disordered states, we think this fact does not constitute a serious drawback. 



C. Topological features of double well potentials 

Double well potentials can be characterized by many quantities: the first one is the 
asymmetry A, that is the energy difference between the two connected minima of a DWP. 
This parameter is essential for the identification of candidate two level systems (TLS), which 
require a value of A of the order of less than f«l K (however, as we will see, this condition is 
not sufficient to identify a TLS). Our results indicate that the asymmetries are distributed 
exponentially, the most part being lower than 5000 K; the distribution is not sensitive to 
the system size. Equally important are the energy barriers V, i.e. the energy differences 
between a minimum and the corresponding saddle. Of course for every DWP there are two 
values of V, one for the relaxation process and one for the activation. The distribution of 
the relaxation barriers similar in shape to that of the A's, but on average the values of V are 
smaller, being significantly present only up to 1000 K. There seems to be little correlation 
between the asymmetry and the barrier (both for the activation and relaxation). 

We have evaluated also the (mass weighted) euclidean distance between pairs of minima, 
defined as 



dist{a, h) 



'yZ~^\^i,a, — '^i,b\'^ , (2) 

, m 



1/2 



where a and b are the two minima, the single atom mass, ri^a its position in the a 
configuration, and m = ^^mi/N. In Fig. 3 we report the relative distributions. The 
distributions do not extend very much beyond ^lOcr and present a maximum at 2a. 
Another very interesting parameter to consider is the participation number defined as 

^part = E#- (3) 



max 



Here di refers to the atom-atom distances in (y) and d^ax is the distance of the atom that 
moves most during the transition; the partecipation number gives an idea on the numbers 
of atoms involved in a transition. In Fig. 4(a) we report this quantity for the three cases 
studied; the same quantity normalized to the number of particles constituting the clusters is 
reported in Fig. 4(b). We see that A^part is a nearly scaling quantity with A^, indicating that 
at least an appreciable part of the atoms that move in the transitions belong to the bulk. 

We can make use of multidimensional transition state theory (for a review see Hanggi 
1985) to establish a link between the potential energy landscape and the (classical) relaxation 
dynamics of the system. Under appropriate conditions (Hanggi 1985) the classical, thermally 
activated transition probability between two metastable minima is given by 
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(4) 



(5) 



where v* is an effective frequency that takes into account the effect of all the degrees of 
freedom on the one which forces the transition. This frequency can be rewritten as 



with i/q the lowest frequency of the dynamical matrix in the starting minimum, while R is 
the product of all the N — 1 positive eigenvalues of the dynamical matrix at the saddle point, 
divided by the corresponding product at the minimum. In Fig. 5 wc show the value of this 
entropic factor R (which enhances or reduces the transition probability), versus the energy 
barrier value, both for the relaxation and the activation transitions. This plot evidences a 
rather marked correlation between R and the barrier height, and in most cases i? <1 or even 



We describe now the most difficult part of our work, that is the selection of two-level 
systems. Following the energy landscape paradigm, these are DWP that imply purely quan- 
tum mechanical relaxation processes. The calculation has been carried out by assuming the 
validity of the ID Wcntzel-Kramers-Brillouin (WKB) approximation (Froman and Froman 
1965, Landau and Lifchitz 1967; Schiff 1968; a discussion on the validity of this approxima- 
tion in the case of Ar clusters, as well as on the effects to be expected when it is released, 
can be found in Daldoss et al. 1999): the splitting of the ground state is given in terms 
of the action integral between the two minima. In principle, in order to apply the WKB 
procedure it would be necessary to find the least action path, i.e. the classical path that 
takes from one minimum to the other and minimizes the action integral; this involves rather 
heavy numerical calculations. As a starting point, we decided to use the path that takes 
from one minimum to the other, and that is defined at each point by the direction of the 
minimum eigenvector, as described in Section II. Work on the minimization of the action 
integral is in progress (Brangian et al. 1999). It should be stressed that as cumbersome as 
this procedure may look, it is probably the simplest way to get a quantitative estimate of 
the tunneling splitting. It should also be mentioned that, in agreement with previous works 
(Heuer and Silbey 1993, Heuer 1997, Daldoss et al. 1998), we find roughly one TLS for 1000 
DWP's, which imphes a very extensive search strategy. 

We have performed an a posteriori test on the reliability of the use of the ID WKB 
scheme. The main hypothesis is that, along the chosen ID path (the least-action or a 
very close one), the degrees of freedom other than the considered one are decoupled from 
it: in this case it is reasonable to expect that the Schroedinger equation may be (nearly) 
factorized into 3N mutually independent equations, a condition for the applicability of the 
WKB scheme in many dimensions (Schiff 1968). If this factorization should actually take 
place, the eigenvalues of the dynamical matrix (which are proportional to the curvatures of 
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D. Two-level systems (TLS) 
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the hypcrsurface) other than the lowest one should remain constant along the chosen path. 
The departure of these eigenvalues from constancy gives a measure (though a qualitative 
one) of the invalidity of the ID scheme. In Fig. 6 we show, as an example, the results for 
two TLS belonging to clusters of 90 atoms. In one case (Fig. 6(b)) the approximation is 
satisfied in a very good way; on the contrary, in the other case (Fig. 6(a)) there is appreciable 
mixing of the low energy eigenvectors. These results are in qualitative agreement with those 
relative to Ar clusters (Daldoss et al, 1999); they imply that the actual structure of TLS in 
disordered systems is probably much more complex than expected, since many-dimensional 
effects seem to play important roles. 

IV. CONCLUSIONS 

We have reported preliminary results of an extensive investigation on the properties of 
the potential energy landscape in Si02 clusters of three different sizes (60,90, and 150 atoms); 
the aim of this research is to find connections between the properties of the landscape and the 
high- and low-temperature relaxation dynamics. By analyzing the structure of the clusters 
and the topological features of their energy landscape (and in particular of its minima and 
first-order saddle points), we have identified tunneling centres and studied the conditions of 
validity of the WKB approximation, which allows a quantitative estimate of the tunneling 
splitting. 
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TABLES 



TABLE I. Comparison between experimental bond lenght obtained in vitreous silica and those 
calculated from our clusters. 



C^cxp [A] d^mni [A] 


St 


-0 


1.61 1.62 





-0 


2.63 2.61 


Si 


-Si 


3.13 3.04 



FIGURE CAPTIONS 



Fig. 1. Radial distribution function g{r) in clusters of three different sizes. Besides the 
three main peaks, which arc also found in bulk systems, we notice other peaks for the 0-0 
and Si-Si bonds, probably due to the presence of surface co-ordination defects. 

Fig. 2. Energy distribution for saddles, upper and lower minima in clusters of 150 atoms. 

Fig. 3. Distribution of euclidean distances between connected minima, in units of A and 
a — I.GA, i.e. the Si-0 bond length. 

Fig. 4. Participation number: the y axis units are chosen such that the integrals of the 
various curves is equal to 1. The bottom plot refers to the participation ratio normalized to 
the number of atoms. 

Fig. 5. Entropic ratio as a function of the barrier height for N =60, 90, and 150. 

Fig. 6. Variation of the 10 lowest eigenvalues of the dynamical matrix along the minimum 
eigenvalue path (see text) in two TLS of clusters of 90 atoms. The minima are arbitrarily 
assigned the ±1 values of the coordinate along the path. 
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Brangian et al, Fig. 2 



10 







distances (a) 
10 15 



480 
320 
160 




o 240 
o 

3 120 

o 

o 




200 

150 

100 h 

50 





1 1 r 



20 



[SiOJ 



2-" 20 




[SiOJ 



2-" 30 



[SiOJ 



2-" 50 



jyyyyyyyyByBiiin.- — 







10 20 30 

o 

distances (A) 



25 




40 



Brangian et al. Fig. 3 



11 



0.10 

0.08 

0.06 

0.04 

0.02 

0.00 
0.06 

0.04 

0.02 

0.00 







10 



participation ratio 
20 30 40 



50 



(a) 



[SiOJ 



□ □ □ 

/ h 

□ \ 



2 J 20 

[SiOJ3, 
-V- [SiOJ^, 



^ I 1 1 1 'ta jT-rti 



1^ 



(b) 



60 




0.2 0.3 0.4 0.5 

normalized participation number 



Brangian et al. Fig. 4 



12 




Brangian et al. Fig. 5 
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